Measuring nonadiabaticity of molecular quantum dynamics with quantum fidelity and 

with its efficient semiclassical approximation 
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We propose to measure nonadiabaticity of molecular quantum dynamics rigorously with the quan- 
tum fidelity between the Born-Oppenheimer and fully nonadiabatic dynamics. It is shown that this 
measure of nonadiabaticity applies in situations where other criteria, such as the energy gap crite- 
rion or the extent of population transfer, fail. We further propose to estimate this quantum fidelity 
efficiently with a generalization of the dephasing representation to multiple surfaces. Two variants 
of the multiple-surface dephasing representation (MSDR) are introduced, in which the nuclei are 
propagated either with the fewest-switches surface hopping (FSSH) or with the locally mean field 
dynamics (LMFD). The LMFD can be interpreted as the Ehrenfest dynamics of an ensemble of 
nuclear trajectories, and has been used previously in the nonadiabatic semiclassical initial value 
representation. In addition to propagating an ensemble of classical trajectories, the MSDR requires 
evaluating nonadiabatic couplings and solving the Schrodinger (or more generally, the quantum 
Liouville-von Neumann) equation for a single discrete degree of freedom. The MSDR can be also 
used to measure the importance of other terms present in the molecular Hamiltonian, such as di- 
abatic couplings, spin-orbit couplings, or couplings to external fields, and to evaluate the accuracy 
of quantum dynamics with an approximate nonadiabatic Hamiltonian. The method is tested on 
three model problems introduced by Tully, on a two-surface model of dissociation of Nal, and a 
three-surface model including spin-orbit interactions. An example is presented that demonstrates 
the importance of often-neglected second-order nonadiabatic couplings. 



I. INTRODUCTION 

Nonadiabatic effects give rise to a great variety of phe- 
nomena in chemical dynamics i^r— To account for these 
effects, many theoretical methods have been developed. 
The most accurate but also the most computationally 
demanding are wave packet approaches which solve the 
Schrodinger equation for both electrons and nuclei di- 
rectly. Some wave packet methods, e.g., the multi- 
configuration time-dependent Hartree method,"*'^- have 
been successfully applied to problems with tens of degrees 
of freedom. Trajectory based nonadiabatic Bohmian 
dynamics^ is another in principle exact method, which 
can be, moreover, combined with electronic structure 
computed on the fly,^ Less accurate but also less expen- 
sive are various semiclassical approaches, which can also 
describe some quantum effects, especially on the nuclear 
motion. These include multiple-spawning methodsj^iii^ 
methods based on the Herman-Kluk propagator f^ir— or 
the surface hopping and jumping method of Heller et al^ 
The most widely used are methods in which the nuclei 
are treated classically and the quantum effects enter only 
through interaction with electrons, which are described 
quantum mechanically. Among these belong methods 
based directly on the mixed quantum-classical Liouville 
equatioufi^"— the mean field Ehrenfest dynamics, vari- 
ous surface hopping methods)^"— or methods in which 
the classical limit is obtained by linearizing the path in- 
tegral representation of the quantum propagator f2§, 

Unfortunately, all of these methods are significantly 
more computationally demanding than their adiabatic 
or diabatic counterparts. (In the following, we discuss 
mostly nonadiabatic dynamics. Nevertheless, the discus- 



sion holds almost entirely also for the nondiabatic dy- 
namics in the diabatic basis.) One goal of this paper is 
to find a general criterion which would determine when a 
given dynamics is nonadiabatic enough to justify the use 
of the expensive nonadiabatic methods. Several possible 
criteria could be envisaged. One widely used criterion is 
the extent of population transfer or, more precisely, the 
decay of survival probability Pqm on the initially popu- 
lated potential energy surface (PES) due to nonadiabatic 
couplings. Although a fast decay of -Pqm is a clear sign 
of nonadiabaticity of the dynamics, the opposite impli- 
cation is not necessarily true, as can be in seen in Ref. 
|29| and as will also be demonstrated below. Another cri- 
terion, often employed to decide when to "switch on" the 
couplings in nonadiabatic calculations, ~ is the energy 
gap criterion: one simply monitors the energy difference 
between PESs and when it becomes sufficiently small, 
the dynamics is considered nonadiabatic. While useful 
in practical calculations, this criterion does not always 
reflect the actual extent of nonadiabaticity, which also 
depends on nonadiabatic couplings and on the nuclear 
momentum. Other approximate criteria, which are inter- 
mediate between the two criteria mentioned, estimate the 
change of Pqm from basic properties of the PESs, from 
couplings between the PESs, and from the nuclear veloc- 
ity. Examples include variants^'^ of the Landau-Zener- 
Stiickelbergnrodel.^^"^ 

In Ref. |29| we proposed a more rigorous quantitative 
criterion of the non(a)diabaticity of the quantum dynam- 
ics, based on the quantum fidelity Pqm'^^ between the 
adiabatically and nonadiabatically propagated molecular 
wave functions. More precisely. 



FQM{t)^\fQM{t)?^\{i'\t)\rm 



(1) 



2 



where is the quantum state of the molecule 

evolved using the adiabatic Hamiltonian H'^ with uncou- 
pled PESs and \^'^{t)) is the quantum state evolved using 
the fully coupled nonadiabatic Hamiltonian H'^. When 
Fqm ~ 1, li^^i^)) is close to \'ip'^{t)) and an adiabatic 
simulation is a good approximation to the nonadiabatic 
simulation. When Fqm <C 1, adiabatic treatment is in- 
adequate and a nonadiabatic method should be used. 
Unlike the energy gap and population transfer criteria, 
the fidelity criterion can detect more subtle nonadiabatic 
effects caused, e.g., by the displacement and/or inter- 
ference on a single PES surface (see Fig.[T|). Panels (c) 
and (d) of Fig.[T] show two extreme examples in which 
the nonadiabatic couplings induce a hop to the excited 
surface followed by a hop back to the ground surface so 
that at large times, only the ground state is occupied. 
While the nonadiabatic couplings have no effect on the 
survival probability on the ground surface, they have a 
large effect on the molecular wavefunction since the re- 
turning wavepacket may accumulate a time delay (panel 
c) or a phase (panel d), and hence can have a zero overlap 
with the wavepacket propagated adiabatically. Although 
neither case can be detected by the survival probability 
criterion, both scenarios can be detected easily by fidelity 
©• 

Due to the generality of definition ([l}, fidelity can be 
employed in many applications with nonadiabatic dy- 
namics, depending on the choice of and iJ*^. Very re- 
cently, this fidelity criterion of nonadiabaticity was used 
to find the optimal time-dependent Hamiltonian maxi- 
mizing the adiabaticity of the dynamics from an initial 
state to a desired target state. '^^ Besides the use of fidelity 
to measure nonadiabaticity, two other applications are 
explored here. First, we consider the importance of addi- 
tional terms in a nonadiabatic Hamiltonian, such as spin- 
orbit coupling terms or couplings to an external field. In 
this case, both if" and H'^ are coupled by nonadiabatic 
coupling terms, and the additional term of interest, miss- 
ing in H^, is present in H"^. Second, Fqm may be used to 
evaluate quantitatively the accuracy of the quantum dy- 
namics with an approximate or interpolated nonadiabatic 
Hamiltonian H'^ in comparison to the quantum dynam- 
ics with an accurate nonadiabatic Hamiltonian iJ*^. This 
application is a generalization to nonadiabatic dynamics 
of the idea proposed for adiabatic dynamics in Refs. [sol 
andiM 

The remaining (but difficult) question is how to com- 
pute -Fqm- The most straightforward way would be to 
propagate the wave functions with some wave packet 
method and to use Eq. ([1]) directly. This approach, how- 
ever, suffers from the previously mentioned disadvantages 
of wave packet methods. Instead, below we derive an 
accurate, yet efficient semiclassical method capable of 
computing not only fidelity -Fqm but also fidelity am- 
plitude /qm- The method, which we refer to as multiple- 
surface dephasing representation (MSDR), is a general- 
ization of the dephasing representation (DR);^"— de- 
rived for the adiabatic dynamics using the Van Vleck 



propagator. In the single surface setting, the DR is 
closely related to the semiclassical perturbation approx- 
imation of Hubbard and Miller— and phase averaging 
of Mukamel.^^ Its main applications include calculations 
of electronic spectra^"— and evaluations of stability of 
quantum dynamics The main advantage of the 
MSDR compared to wave packet methods is that the 
computational cost of MSDR, similarly to the DR. does 
not scale exponentially with the number of degrees of 
freedom!^ The MSDR can therefore be applied to prob- 
lems with dimensionality far beyond the scope of cur- 
rent methods of quantum dynamics. The advantage of 
MSDR in comparison to most other semiclassical ap- 
proaches is that the MSDR does not require the Hes- 
sian of the potential energy, which is often the most ex- 
pensive part of semiclassical calculations (see, e.g., Ref. 
ISSl ). Finally, in contrast to methods treating the mo- 
tion of nuclei completely classically, the MSDR includes 
some nuclear quantum effects approximately via the in- 
terference between the classical trajectories representing 
a wave packet. 



The MSDR is not the first extension of the DR to 
nonadiabatic dynamics. In Ref. [2^ we have intro- 
duced another extension of the DR, which is here referred 
to as integral multiple- surface dephasing representation 
(IMSDR) and which performs satisfactorily in the case 
of nearly diabatic dynamics in the diabatic basis. Unfor- 
tunately, the accuracy of the IMSDR deteriorates when 
the dynamics is far from the diabatic limit. Another im- 
portant limitation of the IMSDR is that it cannot be 
used in the adiabatic basis. Below we shall demonstrate 
that the MSDR is both more accurate and more gen- 
eral than the IMSDR. A small price to pay for this im- 
provement is that in contrast to the IMSDR, in which 
fidelity is computed as an interference integral at the end 
of dynamics, in the MSDR the Liouville-Von Neumann 
equation for one discrete (collective electronic) degree of 
freedom has to be solved during the dynamics. For pure 
states, this equation is simple and is equivalent to the 
Schrodinger equation for one discrete degree of freedom 
which is already solved during the Ehrenfest or fewest 
switches surface hopping (FSSH) dynamics. Both MSDR 
and IMSDR reduce to the original DR for systems with 
a single PES. 



The outline of the paper is as follows: in Sectionllll 
the MSDR is derived. In Section HIIl the method is used 
to evaluate nondiabaticity and nonadiabaticity of quan- 
tum dynamics in model cases, to evaluate the importance 
of an additional coupling term in a nonadiabatic Hamil- 
tonian and to evaluate the accuracy of an approximate 
Hamiltonian. Computational details are summarized in 
the same section. Section lTVl concludes the paper. 
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Figure 1. Possible criteria of non(a)diabaticity of quantum dynamics. The static energy-gap criterion does not take into account 
the dynamics of the wave packet (a). The population transfer criterion measures the actual decay of probability density on 
the initial PES (b). It is more sensitive than the static energy-gap criterion. Fidelity criterion can capture the population 
transfer (b) as well as other nonadiabatic effects such as displacement (c) or interference (d) on a single PES, which would be 
undetected by the population transfer criterion. is the decoupled Born-Oppenheimer Hamiltonian, whereas H'^ is the fully 
coupled nonadiabatic Hamiltonian. 



II. THEORY 



MSDR 



A starting point for the derivation of the MSDR is the 
expression for quantum fidelity amplitude formulated in 
terms of the density matrix^ 

/qm it) = Tr (e-'H'*/'' • p'"'' ■ e+'^^"'/^) , (2) 

where p""* is the density operator of the initial state, 
and are two different molecular Hamiltonians ex- 
pressed either in the diabatic or adiabatic basis. (Bold 
face denotes n x n matrices acting on the Hilbert space 
spanned by n electronic states, hat'denotes nuclear oper- 
ators.) Note that Eq. ([2]) applies to both pure and mixed 
states. Formally, can be written as = + eV, 
where e controls the extent of the perturbation. Express- 
ing /qm in the interaction picture gives 



/qm (t) = Tr p'-'.E(t) 



(3) 



where 



is the echo operato r^^i^^ and 



(5) 



is the perturbation in the interaction picture. A par- 
tial Wigner transform^ of Eq. ^ over nuclear degrees 
of freedom yields an alternative exact expression for the 
fidelity amplitude, 



/qm (t) = h-'^Tr, J dXp^'' {X) ■ Ew {X, t) , 



(6) 



where Aw(^) is the partial Wigner transform of opera- 
tor A, 



Aw (A) = / 



g + lNexpl'ji^), (7) 



X denotes the point (Q, P) in the 2_D-dimensional phase 
space, and Trg is the trace over electronic degrees of free- 
dom. Direct evaluation of the Wigner transform of the 
echo operator is unfortunately difficult without approx- 
imations. As the first approximation, one may truncate 
the Taylor expansion of the exponential operator in the 
Wigner transform of a product of two operators 



(A.B) 



w 



exp 




= Awexp B^^8) 



after the zeroth-order term. In Eq. ([5]), {A, B} 



^ • ^ — ^ • ^ is the Poisson bracket over the nu- 
clear degrees of freedom and arrows indicate on which 
argument the derivatives act. An iterative application 
of expansion (|8]) to the echo operator (|4]) expressed as a 
time-ordered product of the short time propagators gives 

To evaluate this expression, the time evolution of 
(A, t) is required. The second approximation in- 
volves replacing the exact equation of motion by a mixed 
quantum-classical (MQC) propagation scheme described 
below, leading to the final expression for MSDR of fidelity 
amplitude, 

/msdrW = h-^'TTeJ dXp^^\X) ■ Te-'-^o^'w.MQcix,t')dt'/n 
where the average in the last row is defined in general as 



(A (A)) 



Tre / dXpjX) ■ A (A) 
Tre / dXp{X) ■ 



Equation (jlO[) makes it clear that the only fundamental 
difference between the MSDR and IMSDR introduced in 
Ref.,23 is the time ordering operator T present in the 
MSDR but missing in the IMSDR. 
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B. Propagation scheme 



Equation (fTO| gives /] 



/MSDR 



in terms of VJ 



W,MQC 



what remains to be done is finding a prescription 
for ]y[QQ(X, t). The exact equation of motion for 
{X, t) in the interaction picture is 



dV\v iX,t) 
dt 



H°,Vi(t) {X) 



w 



(11) 



Note that Eq. pip is, up to the sign, the same as the par- 
tiaUy Wigner transformed LiouviUe-Von Neumann equa- 
tion describing the propagation of the density matrix of 
the unperturbed system: 



dp.^{X,t) 
dt 



w 



(12) 



We will take advantage of this analogy and derive our 
approximate propagation scheme from Eq. (|12p instead 
of from Eq. (fTT|) . This will turn out to be slightly easier 
and, more importantly, we will simultaneously find an 
approximate solution of a much more general problem of 
propagating the density matrix. Below we omit super- 
script in H*^ since now we deal with a single Hamil- 
tonian. In the system consisting of light electrons and 
heavy nuclei, Eq. (IT^ is approximated to the first order 
in the mass ratio by the mixed quantum-classical 

Liouville equatioii* 

^Pw,MQC 



M 
17.58-63 



dt 



— ^ T [Hw,Pw,mqc] 



+ 2 ({Hw,Pw,MQc} ^ {Pw.MQC'Hw}) , 

(13) 

where the explicit dependence of mqc '^^ time and 
on the nuclear phase space coordinate X was omitted for 
clarity. 

Equation ([TU)) together with Eq. define the 

MSDR. Several numerical approaches exist that solve 
Eq. (fT3|) in terms of "classical" trajectories X{t). How- 
ever, since trajectory based methods for solving Eq. 
are still relatively complicated, below we study two vari- 
ants of MSDR where Eq. ([T^ is further approximated. 
The common feature of the two approximations is that 
all elements of pyj^ (X, t) are propagated using the same 
PES (which may, nevertheless, differ for different trajec- 
tories). For simplicity, from now on the subscript MQC 
is omitted. 



1. Locally mean field dynamics 

The first approach starts by rewriting [X, t) as 

p^{X,t) ^ p{X,t)p,{X,t), (14) 

where p{X,t) := Trep^(X, t) is a scalar function of X 
and t and hence TrePg(X, t) — 1 for all X and t. By 



substituting the still exact ansatz into Eq. (IT511 . one 
obtains 



dt^^ 



' at 



-p[Hw,Pe] 



ldp_ 
2dP 
1 
2 

ldp_ 
2dQ 



w 



dQ 

dQ '~dP 
9Hw 



dP 



w 



dP 



dQ 



(15) 



where [A, B]^ = A • B -I- B • A is the anticommutator. 
In the next step, the trace over the electronic degrees of 
freedom is performed, which in the diabatic basis leads 
to the following equation of motion for p {X, t) : 



dp _ dp 1 9Hw \ dpi 9Hw 
dt^dP\ dQ /^~dQ\ dP 



V dQ 



(16) 



where (A)^ = Trg (p^ • A) is a partial average of A over 
the electronic subspace and where we have used that 



Tr, 



9Hv 



' dQ ) 



rTr, 



0. (The equation 



^ dP dQ J ~ M 
of motion in the adiabatic basis will be derived later.) 
Substitution of Eq. ([TE)) back into Eq. (IT5|) yields 



dPe « rrr 









dQ 




an 



w 



dQ 





'dUw 




dP 



an 



w 



dP 



where identity ^ 



dP ' dQ 



P_9p, 
M dQ 



/9Hw\ 
\ dP U 



(17) 



3Pe 

dQ 



was used in the fifth row. Both Eqs. (|16l) and pT)) contain 
terms of the form {^^)^§^ and ^ 



dQ 



1^, which 



can be combined with the time derivative ^ to form the 
convective derivative 

£l^^ + d^ + p^ 

Dt dt ^ dQ dP 



5 



and transform the equations from the Eulerian reference 
frame at rest to the Lagrangian reference frame moving 
with the phase space flow given by the vector field 



Q,P] = 



9Hv 



9P \ dQ 

In the Lagrangian frame, Eq. (|16p transforms to 



Dt ^ " y dQ ' dP 



(18) 



(19) 



Since the two terms in the fourth row of Eq. (IT71) cancel 
each other exactly, this equation transforms to 



Dt 







1 


[ dQ '^"^J 





'dUw 


dPe 




_ dQ 


dP 




PeTre 



(9H 



w 



dQ 



dPe] 
dP I 



(20) 



Note that the second and third row of the right hand 
side of Eq. ((20)) contain differences (put in parentheses 
for emphasis) between products of the electronic density 
matrix (or its P derivative) with the nonaveraged and 
averaged gradients of the Hamiltonian. Therefore, the 
second and third tow may often be small compared to 
the first and fourth row of Eq. (|20l) . The last term on 
the right hand side of Eq. (piU) corresponds to the right 
hand side of Eq. (jl9p and is responsible for maintaining 
TrePe(X, t) = 1 during the (non-approximated) MQC 
dynamics. Until now all operations have been exact and 
Eqs. and ([^0)1 are equivalent to the original mixed 
quantum classical Liouville equation p3p . 

Now we will make the locally mean field approxima- 
tion: Specifically, once in the Lagrangian reference frame, 
all terms in both Eqs. and (PU)1 containing phase 
space derivatives of p (X, t) or (X, t) are neglected to 
obtain the approximate locally mean-field equations of 
motion. The resulting equation for p {X, t) in the La- 
grangian reference frame is simply 



Dp 



LMFD 



Dt 



0, 



(21) 



and the equation of motion for the electronic part of the 
density matrix {X, t) is 



Dp, 



e,LMFD 



Dt 



- [Hw, Pe,LMFDj 



(22) 



Both equations can be combined together and trans- 
formed back to the Eulerian reference frame to yield the 



equation of motion for the total density matrix p^(X^ t), 

,LMFD 



dt 



Hw, Pw.lmfdJ 
9Hw 



h 

.LMFD 

dP \ dQ 

^Pw.LMFD P 



dQ 



M' 
_p 

M- 



(23) 



where we have used that ( ) 

We call the dynamics expressed by Eq. (|23p [or, equiv- 
alently, by Eqs. pB|) . (PT|) and ([^ j the locally mean field 
dynamics (LMFD) since the force acting on nuclei at po- 
sition X is the force averaged over the "electronic" part 
of the density matrix p^{X,t) at X. Note that the well- 
known Ehrenfest dynamics (which uses a single nuclear 
trajectory) can be derived in a similar way using an ap- 
proximate ansatz = 5{X — X {t))p^ (t), where S is 
the Dirac delta distribution and X (t) is the phase space 
coordinate at time t.-^ Similarly, a truly mean field dy- 
namics for a wave packet different from a 6 distribution 
can be derived using an (again approximate) ansatz in the 
form of a Hartree product p^™ {X, t) = p{X, t)p^ (t). 

As can be seen from Eq. (P^ . for pure states, each 
phase space point is propagated by the mean field Ehren- 
fest dynamics, according to Eq. ([TS)) . Nevertheless, 
in contrast to the purely mean field dynamics, Eq. (|23p 
describes the propagation of the density p^{X,t) us- 

and (^1^)^ for each 



different values of 



mg ^..xv..^... \ -gQ-/ \^p- 

value of X . Interestingly, this dynamics corresponds ex- 
actly to the nuclear dynamics appearing in the nonadia- 
batic IVR model^iS^ which uses the Meyer-Miller-Stock- 
Thoss Hamiltoniani^iS Also note that outside of cou- 
pling regions and when only one surface is occupied, re- 
sulting trajectories are equivalent to those obtained with 
the Born-Oppenheimer dynamics. 

To derive a corresponding LMFD equation of motion in 
the adiabatic basis, one can express the mixed quantum 
classical Liouville equation in the adiabatic basis,— 
use the exact ansatz , and apply a similar LMFD ap- 
proximation as above. However, this procedure is quite 
tedious in the adiabatic basis. The LMFD in adiabatic 
basis can be obtained more easily by directly transform- 
ing the final LMFD equation of motion (PT|) from diabatic 
to adiabatic basis using the relation&ii 



dA 



w 



dQ 
dP 



^-[A^,F] and 
dP ' 



(24) 
(25) 



where the superscript A denotes the transformation to 
the adiabatic basis and F is the vector matrix of nonadia- 
batic coupling vectors. Specifically, A^(X, t) is a matrix 
obtained by first partially Wigner transforming a general 
operator A(t) [to form Aw(A", t)] and then by trans- 
forming Aw (AT, t) into the adiabatic basis. Matrix F is 
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defined componentwise by Fjk (Q) = {ctj {Q) I^Q^fc iQ))-. 
where \ak {Q)) is fc-th element of the adiabatic basis at 
position Q. Applying relations and (P5)l to deriva- 
tives of both Hw and lmfd Eq. ([23]) immediately 
yields the LMFD equation of motion in the adiabatic ba- 



trajectory according to the following prescription: For 
a trajectory starting at X, the probability for its initial 
surface to be surface j is given by the diagonal element 
Equation pn]) for fidelity amplitude is rewrit- 
ten as 



SIS. 



dpi ,LMFD 

dt 



/msdr (i) 



W-«^^F(Q),PW,LMFD 



{X) 



(28) 



dpi ,LMFD 
dP 

^Pw.LMFD P 



dnl 

dQ 



— ( [H^, F] ) j where we have used the notation 



dQ 



M' 



(26) 



where is the diabatic Hamiltonian matrix Hw ex- 
pressed in the adiabatic basis. Note that consists 
only of the kinetic energy term and the diagonal adia- 
batic potential energy matrix; in particular, does 
not contain the nonadiabatic couplings. Again, the dy- 
namics of a single trajectory is identical to the Ehrenfest 
dynamics. 



_ JdXp{X)A{X) 
(^(^))p(x) - JdXpiX) 

and the fact that (Xrepj."" (^))p.nH(x) = 1- The time- 
ordered product T e'^-Zo' '^w('^.*')''*'/'i jg evaluated along 
each trajectory by successive multiplication of short time 
propagators corresponding to each time step. The expo- 
nent of each short time propagator is computed using 

YliX, t) = UO(X, t)-i . Vw(XO(i)) • UO(X, t), (29) 



2. Fewest switcties surface hopping 

The second alternative approximate propagation 
scheme is based on the physically motivated FSSH 
algorithm. In this scheme, each phase space point rep- 
resenting the Wigner density distribution is propagated 
independently using the FSSH dynamics. Compared to 
the LMFD, the FSSH is used at no additional cost, ex- 
cept for the stochastic hopping algorithm itself, because 
the same Eq. ([^^ (or its equivalent in the adiabatic ba- 
sis) has to be solved during the electronic part of the 
dynamics. On the other hand, averaging the force over 
electronic states is avoided in the FSSH. As will be shown 
below, neither method is universally optimal and the best 
propagation method depends on a problem studied. 



C. Algorithm 



1. General initial state 



where 



(30) 



and X^{t) denotes a trajectory of starting at 

A'°(0) — X . Equation can be vaguely interpreted as 
a combination of a quantum interaction picture for the 
electrons and a "classical interaction picture" (in which 
the perturbation is neglected) for the nuclei. Operator 
U° {X, t) is again computed by successive multiplication 
of short time propagators along the trajectory. At the 
end of the dynamics, the weighted phase space average 
in Eq. is computed as the arithmetic average over all 
trajectories: 



/msdr {t) 



N=l 



pr'(^iv)-r n 



-i£V^(Xjv,MAi)At/'i 



M=l 



(31) 



To compute /MSDR(i), we rewrite the initial density 
matrix exactly in the form 

p^"(X) = p™*(^)pf*W, (27) 

where p'"'' {X) :— TrePw' (^)- Scalar nuclear density 
/9'"'*(X) is sampled by iVtraj phase space points that serve 
as initial conditions for A'traj trajectories propagated us- 
ing either the LMFD or FSSH dynamics. For each gener- 
ated phase space point, the electronic part p™'* {X) satis- 
fies Trepj."" {X) = 1. In the case of LMFD, X determines 
the initial condition completely. In the case of FSSH, one 
also needs to select the initial surface randomly for each 



2. Initial Hartree product state 

The algorithm simplifies slightly when the initial den- 
sity operator p'"'* is a Hartree product 

p""* = /5i^;'*(»pr*, (32) 

where /5™'* and p™'* are the nuclear and electronic den- 
sity operators, respectively. This Hartree approximation 
is usually an excellent approximation outside of coupling 
regions when only one PES is occupied. For initial den- 
sity in the product form dSl]), p™* {X) = p'"" (X) p|f". 
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where p'"'* (X) =JPw']w(^) P^'^ independent 
of X. Equations (P8|) and pTjl become 



/i 



MSDR 



(33) 



and 

/msdr (<) 



_init 



iVtr 



-ieV^(XN,AfAt)At/?i 



JV=1 M=l 



(34) 



where (X, t) is the wave function with initial condition 
{X, 0) = c'"'* (X) that solves the Schrodinger equation 
^2^1^ -iH^(XO(i)) • c<^ (X,t) for a single discrete 
electronic degree of freedom in the Lagrangian reference 
frame given by H^. In both propagation algorithms cur- 
rently used with the MSDR, i.e., the LMFD and FSSH 
dynamics, c° {X, t) is already available; only c*^ (X, t) has 
to be determined additionally. Expressed explicitly in 
terms of trajectories, Eq. ([39l) allows computing the fi- 
delity amplitude simply as 



N, 



/MSDR 



it) 



traj 

\ " „0 

N 



^traj 



(X,t)^c5v(X,t). (40) 



3. "Electronically pure" initial state 



III. RESULTS 



The algorithm simplifies also for "electronically pure" 
initial states, by which we mean (in the most general 
sense) states for which the electronic density matrix 
pimt i-jjg product (|?7)) is pure for all X, i.e., satis- 



fies the condition Trg 
as the tensor product 



{xy 



1 and can be written 



(x)®c™* {xy 



(35) 



in terms of an initial electronic wave function c""* (X). 
The generalized electronically pure states include, as a 
special case, the Hartree product (|32|) in which the con- 
stant electronic matrix Pg"'* is pure (while the nuclear 
factor p™' may be mixed). To derive the simplified ex- 
pression for /msdr (t) , we first rewrite the approximate 
Wigner transform of the echo operator in Eq. (|28| as a 
product of the electronic evolution operators: 

^^~^eJ,';^r'„{x,t')dt'/h ^ u" {X,ty^ ■ U'{X,t), (36) 

where U°(X, t) was defined in Eq. ([50)) and, similarly, 

V%X,t) = 7-e-*/o HW(^''(t'))rft'A, (37) 

Note that the electronic evolution operators {X, t) 
and (X, t) are defined separately for each nuclear tra- 
jectory. Using expression p6l) , we can rewrite the MSDR 
of fidelity amplitude ([25)1 as 



/msdr {t) 



Tr, 



(X)-U" {X,t)~^ -WiX^t) 



(38) 



This expression, which is still valid for any initial state, 
is equivalent to Eq. ^ when either the LMFD or FSSH 
dynamics is used for propagation. 

For electronically pure states Eq. ([55]) further 

simplifies into the weighted phase space average 



/msdr (t) = (c° (X, t)^ -c' {X, t)) ^^^^ 



(X) 



(39) 



To test the MSDR, we considered several model sys- 
tems allowing comparison with the exact quantum- 
mechanical solution in both the diabatic and adia- 
batic bases. Specifically, we used variants of the three 
one-dimensional model potentials introduced by TuUy^^ 
and the simple two-level model of photodissociation of 
jgj^j^68-7o rj]^g j^g^gg -^^ Tully's models was set to 2000 a.u., 
approximately corresponding to the mass of the hydrogen 
atom. The reduced mass in the model of photodissocia- 
tion of Nal was set to 35480.25 a.u. 

The initial density matrix was in all cases a density 
matrix of a pure state, so /msdr was evaluated using 
Eq. (gOl)- In Subsections lIII Al and llll B[ fidelity is used as 
a quantitative measure of the importance of the diabatic 
or nonadiabatic couplings on the dynamics. In other 
words, Hamiltonian H*^ is the diagonal diabatic (Sub- 
sec. 1111^ or adiabatic fSubsec. lIIIBI) Hamiltonian and 
Hamiltonian H*^ is the full nondiabatic (Subsec. lIII A[) or 
nonadiabatic (Subsec. lIII B[) Hamiltonian, respectively. If 
not mentioned otherwise, the dynamics on H° uses the 
LMFD algorithm. Since the PESs of H° are decoupled 
and (with one exception shown in Fig. [3]) only one PES 
is occupied initially, the dynamics reduces to the simpler 
Born-Oppenheimer classical dynamics of an ensemble of 
phase space points. In Subsec FlH C[ more general Hamil- 
tonians and perturbations are studied. 



A. Nondiabaticity of quantum dynamics 

In this Subsection, the fidelity criterion of nondiabatic- 
ity of the quantum dynamics in the diabatic basis is 
studied together with the IMSDR and MSDR approx- 
imations of fidelity. The IMSDR was already shown to 
perform well only when the dynamics was close to the di- 
abatic limit. ?^ On the other hand, as demonstrated here, 
MSDR performs well not only in a broader range of prob- 
lems but also for lower wave packet energies. Compar- 
isons of both methods with numerically exact quantum 
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fidelity (-Fqm) in the diabatic basis can be found in Fig. (3] 
Figure[2](a) shows results for the single avoided crossing 
model of Tully.— The initial wave packet has high kinetic 
energy so that the dynamics is fairly close to the diabatic 
limit. Thus both extensions of the DR work very well; 
^MSDR is almost indistinguishable from Fqm- Very sim- 
ilar results were obtained for the other model potentials 
when the wave packet had sufficient energy to cross the 
coupling region almost diabatically. Figure[2](b) demon- 
strates, on the extended coupling model of TuUy,^^ that 
the survival probability Pqm itself is not always a good 
indicator of nondiabaticity of the dynamics. Here, Fqm 
decays quickly to zero despite Pqm staying close to unity. 
Indeed, the quantum dynamics on and H'^ are very 
different. Both extensions of the DR reproduce the de- 
cay of Fqm very accurately. Using the double avoided 
crossing model of Tully^ Fig.[l](c) shows that MSDR 
can accurately reproduce the fidelity behavior even far 
from the diabatic limit. Not surprisingly, the IMSDR 
method fails here. For comparison. Fig. [21(c) also shows 
two MSDR results obtained by exchanging the roles of 
H" and Since, in contrast to H°, H*^ is coupled, 

both the LMFD and FSSH dynamics allow transitions 
between PESs: both are good approximations of Pqm- 
Finally, Fig.[2](d) demonstrates that because the MSDR 
is a semiclassical method based on classical trajectories, 
not permitting tunneling, the method has to be applied 
with care. In the case shown, far from the diabatic limit, 
the wave packet on H° reflects back from the coupling 
region. In the MSDR this reflection results in "rephasing" 
of the trajectories leading to the rise of fidelity back to 
values close to unity, in disagreement with the quantum 
result. Even if roles of H° and are exchanged, prob- 
lems persist. In the quantum dynamics, a major part of 
the wave packet on the coupled Hamiltonian H'^ passes 
through the coupling region with the aid of diabatic cou- 
plings. When the LMFD is used with MSDR, this be- 
havior and associated fidelity decay are captured quali- 
tatively. When the FSSH dynamics is used, the dynam- 
ics is exactly identical to the dynamics on the uncoupled 
Hamiltonian because all hops in the FSSH algorithm 
are frustrated. This points out the importance of tunnel- 
ing in this relatively narrow region of wave packet ener- 
gies. When the initial kinetic energy is smaller than that 
shown in Fig.[5](d), most of the wave packet bounces back 
even during quantum dynamics on H'^ and the MSDR fi- 
delity approaches the quantum result. When the energy 
is somewhat higher, more trajectories pass the barrier (or 
fewer hops in the FSSH algorithm on H'^ are frustrated) 
and as a consequence the MSDR approaches again the 
quantum result. Note that the relatively good result of 
the IMSDR is accidental since the method is not expected 
to work well so far from the diabatic limit. 



B. Nonadiabaticity of quantum dynamics 

In this Subsection, the fidelity criterion of nonadia- 
baticity of the quantum dynamics in the adiabatic basis 
is studied together with the IMSDR and MSDR approx- 
imations of fldelity. As can be seen in Fig.|31 in the adi- 
abatic basis the IMSDR does not work at all, whereas 
the MSDR works as satisfactorily as in the diabatic ba- 
sis. Figure[2](a) shows results for the photodissociation 
dynamics of Nal close to the adiabatic limit. A low en- 
ergy wave packet oscillates on the excited adiabatic PES, 
crossing the coupling region twice per period, each time 
losing a bit of the excited PES population due to the tran- 
sition to the ground PES. The MSDR reproduces the as- 
sociated fldelity decay very well, whereas the IMSDR fails 
completely. Dynamics even closer to the adiabatic limit 
is shown in Fig.|3](b) using Tully's single avoided crossing 
model A. First of all, note that Pqm is again substantially 
different from the survival probability Pqm- Second, the 
flgure also shows that in certain cases the MSDR based 
on the first order nonadiabatic mixed quantum-classical 
equation [Eq. (|13p ] does not accurately reproduce Pqm- 
Interestingly, the agreement can be improved with the 
second order dynamics, which employs the complete par- 
tially Winger transformed nonadiabatic Hamiltonian 

H^™-P = H^-z;i|:F-|lF-F. (41) 

Strictly speaking such an approach goes beyond the 
MSDR defined by Eqs. ^ and ([13]). Nevertheless, in 
the special case of a decoupled unperturbed Hamiltonian 
H° with only one occupied PES, the second order equa- 
tion of motion can be easily used instead of Eq. ([T3l) , be- 
cause the increased order of dynamics affects only phases. 
Trajectories, which are propagated with the decoupled 
Hamiltonian H'', stay unaffected since the contribution 
from the second order Poisson brackets vanishes. The 
second-order Hamiltonian ()41[) . however, cannot be de- 
rived by generalizing the approach used to derive Eq. 
P6p . In that case, the mixed quantum-classical equation 
()13p was derived from the Liouville-Von Neumann equa- 
tion in the diabatic basis and then transformed into the 
adiabatic basis using Eq. ([M)) to yield Eq. ([^^ . Gener- 
alizing this approach to the second order, no correction 
is obtained to the first order Hamiltonian acting on the 
electronic degree of freedom in the Lagrangian frame. In- 
stead, we have used a generalization to the second order 
of the approach from Ref. fisl . where the mixed quantum- 
classical equation was obtained directly by Wigner trans- 
forming the Liouville-Von Neumann equation expressed 
in the adiabatic basis. The exact reason for the discrep- 
ancy between the two approaches is not yet clear to us 
and will be a subject of further investigation. (Some dis- 
crepancy is actually present already in the first order.) 
In other examples that we have studied, the second or- 
der correction does not significantly influence the results. 
Finally, Fig.[3](b) also demonstrates the convergence of 
the MSDR with the number of trajectories showing that 
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Figure 2. Nondiabaticity of quantum dynamics. Panels (a)-(d) compare the numerically exact quantum fidelity (-Fqm) with 
extensions of the DR in the diabatic basis and with the quantum survival probability (Pqm)- (a) Dynamics close to the diabatic 
limit, (b) Dynamics very far from the diabatic limit, (c) Dynamics in the intermediate range, (d) Dynamics in the region where 
quantum tunneling is important. The insets show the diabatic PESs Vjj{r), diabatic coupling Vi2{r) as well as initial [V'(O)] 
and final [ip"{tf) and il^'^(tf)] wave functions evolved with and H'^, respectively. 



as few as 32 trajectories suffice for an accurate approx- 
imation of i^QM- (The convergence is similar for both 
the first order and the second order dynamics.) For the 
£)j^43,54 g^^-^^ IMSDR,^^ an exact formula exists for the 
number of trajectories as a function of the statistical er- 
ror and of fidelity. For the MSDR, such an exact formula 
has not been derived, but similarly to the DRHi^ and 
IMSDRjS more trajectories have to be used for lower fi- 
delity. Figure[3](c) shows that the MSDR usually retains 
its accuracy even far from the adiabatic limit. Neverthe- 
less, similarly to the diabatic basis, special care should 
be taken in such cases. 

In all examples discussed above, only one PES was oc- 
cupied initially. Yet the MSDR works also with more gen- 
eral initial conditions, as shown in Fig.[3](d) using Tully's 
single avoided crossing model. Here 75 % of the initial 
density is located on the lower PES, and the rest on the 
upper one. Note that when more than one PES of is 
occupied, one must watch out for the intrinsic deficiencies 



of the underlying dynamical methods. If the LMFD is 
used for propagation (not shown), each trajectory prop- 
agates on an average PES, given by the weighted average 
of all occupied PES, even outside of coupling regions. 
When the FSSH based algorithm is used [as shown in 
Fig.[3](d)], other problems may occur in a similar situa- 
tion, i.e., when wave packet dynamics on the lower and 
upper PESs differ substantially outside of coupling re- 
gions. 



C. Importance of an additional interaction in the 
Hamiltonian and accuracy of an approximate 
Hamiltonian 

Now we consider two more general applications of the 
MSDR. In the first application, the MSDR is used to eval- 
uate the importance of an additional interaction term in 
the Hamiltonian. In the second application, the method 
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Figure 3. Nonadiabaticity of quantum dynamics. Panels (a)-(d) compare the numerically exact quantum fidelity (-Fqm) with the 
extensions of the DR in the adiabatic basis and (when applicable) with the quantum survival probability (Pqm)- (a) Dynamics 
close to the adiabatic limit, (b) Convergence of the MSDR and importance of the second order dynamics close to the adiabatic 
limit, (c) Dynamics very far from the adiabatic limit, (d) A more general initial condition with both PESs initially occupied. 
The insets show the adiabatic PESs Vjj(r) as well as initial [V'(O)] and final [il)^{tf) and ip''{tf)] wave functions evolved with 
H'^ and H*^, respectively. 



is employed to evaluate the accuracy of quantum dynam- 
ics on an approximate non(a)diabatic Hamiltonian. In 
both cases, H° and H'^ are coupled and may differ in 
both diagonal and coupling terms. Because even H° is 
coupled, the trajectories used in the MSDR do not any- 
more correspond to simpler Born-Oppenheimer trajecto- 
ries even when only one PES is occupied initially . 

1. Importance of an additional interaction in the 
Hamiltonian 

The importance of an additional interaction in the 
Hamiltonian is tested in Fig.|3](a). This case is simi- 
lar to previously studied problems because and 
differ again only by the presence of a coupling element; 
the difference is that H° is now coupled. Both diabatic 
Hamiltonians contain three PESs, of which the lower two 



are identical to the PESs of Tully's single avoided cross- 
ing model. The third PES is flat with constant energy 
E = 0.15 a. u. In both Hamiltonians, the lower two PESs 
are coupled by the same coupling term V12 as in the orig- 
inal single avoided crossing model. Additionally, in H*^ 
(but not in H"), the highest two PESs are coupled with 

V23 = V*^ - (1 + [C exp i-DQ^)] , (42) 

where C = 0.005 and D = 1.0. A more general complex 
form was chosen to emulate the presence of spin-orbit 
coupling terms, which may also be complex valued. As 
can be seen from the figure, the MSDR reproduces the 
exact decay of fidelity accurately. Note that fidelity de- 
cays significantly even though only approximately 1 % 
of the probability density ends up on the highest PES. 
In other words, this is yet another example, where the 
survival probability is a poor measure of the importance 
of couplings between PESs. 
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2. Accuracy of an approximate Hamiltonian 

The last application of the MSDR that we consider is 
the evaluation of the accuracy of an approximate Hamil- 
tonian. An additional difficulty related to this applica- 
tion is that the electronic basis sets used to represent the 
two Hamiltonians may be different. In the adiabatic ba- 
sis set, electronic basis functions are determined at each 
space point by the Hamiltonian itself. The diabatic ba- 
sis, on the other hand, is usually chosen to diagonalize 
the Hamiltonian at some fixed nuclear configuration. To 
avoid difficulties, one may express both Hamiltonians in 
the same basis set, but this is not always feasible (at 
least not easily). Such problem occur, e.g., when one 
compares two Hamiltonians computed on the fly using 
different electronic structure methods. Since Hamiltoni- 
ans used in electronic structure methods are usually ap- 
proximate, the adiabatic basis functions are not necessar- 
ily identical. Moreover, in some methods, such as DFT, 
electronic basis functions do not even have to be deter- 
mined. Therefore, as we are interested in the dynamics of 
nuclei and in the extent of electronic transitions rather 
than in approximations underlying electronic structure 
methods, we solve this problem by replacing our defini- 
tion of fidelity i^QM with an alternative definition -Fqmj 
in which the overlap matrix of electronic basis functions 
is always assumed to be the identity matrix. 

An "approximate'' Hamiltonian H*^ was created by per- 
turbing one of the parameters in the analytical formula 
for H'' in the diabatic basis. (The perturbation was 
transformed into the adiabatic basis when the dynam- 
ics was done in the adiabatic basis). First, the effect 
of perturbing the slope of PESs in the coupling region 
was studied using TuUy's single avoided crossing model in 
the diabatic basis [see Fig.|3](b)]. The perturbation con- 
sisted in increasing the value of parameter B in Eq. (21) 
in Ref.f25l by 50%. Although electronic transitions are 
significant even on H° itself, the MSDR works very well 
with both the LMFD and FSSH dynamics. Figure|l](c) 
shows the decay of fidelity due to an increased depth of 
the PES well in Tully's double avoided crossing model. 
Calculations were performed both in the diabatic and 
adiabatic basis. [Note that Fq^ may differ between adi- 
abatic and diabatic basis, but only in coupling regions. 
In the case shown in Fig.|3](c) the difference is not sig- 
nificant.] To increase the depth of the well in H'^, value 
of parameter A in Eq. (23) in Ref.i25| was increased by 
10 %. As can be seen in Fig.|3](c), the quantum result is 
reproduced very well by the MSDR based on the LMFD 
in both basis sets. When the FSSH dynamics is used, the 
result depends strongly on the basis. In the adiabatic ba- 
sis, the MSDR reproduces the quantum result quite well, 
whereas in the diabatic basis, the method fails to follow 
-^QM even qualitatively. Worse performance of the FSSH 
dynamics in the diabatic basis is well knowu;^ and in 
this specific case may be attributed mainly to the fact 
that 40% of hops in the diabatic basis are frustrated. 

Figure|3](d) shows an opposite case, where the MSDR 



works better with the FSSH than with LMFD. The model 
used is Tully's extended coupling model in the adia- 
batic basis. To introduce the perturbation, the coupling 
strength in H'^ (parameter B in Eq. (24) in Ref. 25) was 
increased by 50%. The MSDR reproduces the initial de- 
cay of fidelity irrespective of the dynamics used. By con- 
trast, the subsequent rise of fidelity caused by the partial 
reflection of the wave packet could be reproduced only 
with the FSSH dynamics because no reflection occurs in 
the mean field description. 



D. Computational details 

All quantum calculations in the diabatic basis set 
were performed using the second order split-operator 
algorithm. ^•^ Calculations in the adiabatic basis were 
done either by transforming the quantum state into the 
diabatic basis, propagating there, and transforming back 
into the adiabatic basis, or directly with the second 
order Fourier method^ The LMFD or FSSH dynam- 
ics were done using the second order symplectic Verlet 
integrator.^ » 



IV. DISCUSSION AND CONCLUSIONS 

Results presented above demonstrate that quantum 
fidelity is useful as a quantitative measure of nondia- 
baticity or nonadiabaticity of quantum dynamics. More- 
over the MSDR, a semiclassical approximation developed 
in this work, has proven to be a reliable yet very ef- 
ficient substitute for the expensive exact quantum dy- 
namics calculation of fidelity. The MSDR, similarly to 
the less accurate IMSDRj^^ is a generalization of DR 
to non(a)diabatic dynamics. In addition to quantum 
effects originating from the interaction of nuclei with 
electrons, which are included in most mixed quantum- 
classical methods, the MSDR includes also quantum ef- 
fects originating from the interference between mixed 
quantum-classical trajectories representing the evolution 
of the initial density matrix. Two approximate variants 
of the MSDR were developed and studied numerically: 
the former uses the LMFD, derived here, while the latter 
employs the FSSH dynamics. The LMFD, although ob- 
tained in a different way, is equivalent to the nuclear dy- 
namics appearing m the nonadiabatic IVR^i^i^^ which 
is, in turn, nothing else than the Ehrenfest dynamics ap- 
plied separately to each classical phase space point rep- 
resenting the initial density matrix. Both FSSH and the 
Ehrenfest method are relatively simple and often used, 
thus both variants of MSDR can be easily implemented 
into any FSSH or Ehrenfest dynamics code, especially for 
pure states. 

Several applications of the MSDR in nonadiabatic dy- 
namics were presented. First, the method was used to 
approximate a rigorous quantitative measure of nondi- 
abaticity or nonadiabaticity of the dynamics based on 
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Figure 4. (a) Fidelity as a measure of the effects of additional (e.g., spin-orbit) couplings on the quantum dynamics, (b)- 
(d) Fidelity as a measure of the accuracy of quantum dynamics on an approximate Hamiltonian. Comparison of the numerically 
exact quantum fidelity {Fqm or Fqm ) with extensions of the DR. The insets show the diabatic or adiabatic PESs Vjj (r) , diabatic 
couplings Vi2{r) and V23{r) as well as the initial ['(/'(O)] and final [i^^itf) and i^'^{tf)] wave functions evolved with H'^ and H"^, 
respectively. 



quantum fidelity. As such the MSDR may be used to de- 
cide - before running the quantum dynamics itself - which 
PESs or Hamiltonian terms are important. Second, the 
method permits establishing the relative importance of 
several interaction terms in a Hamiltonian. Third, gener- 
alizing one of the applications of the original DRj^Si^ the 
MSDR may be used to evaluate the accuracy of quantum 
dynamics with an approximate non(a) diabatic Hamilto- 
nian. Apart from these applications, the MSDR could be 
used, in principle, to compute all quantities expressible in 
terms of quantum fidelity or quantum fidelity amplitude, 
such as various spectra. 

In Section Um we have demonstrated that for one di- 
mensional model systems the MSDR works often satis- 
factorily even far from the diabatic or adiabatic limit. 
The method has yet to be tested in multi-dimensional 
systems. Nevertheless, results obtained with the origi- 
nal Born-Oppenheimer DR method demonstrate that the 
convergence of the method is actually independent of the 



dimensionality of a problem^ and that the accuracy does 
not deteriorate with dimensionality^- Thus, we expect 
that the MSDR would perform well especially in chaotic 
multi-dimensional systems such as some molecules, pro- 
vided that the underlying mixed classical-quantum dy- 
namics is a reasonable approximation to the quantum 
dynamics. Unfortunately, this is not always the case. 
It is well known that the Ehrenfest dynamics (and also 
the LMFD) can be qualitatively incorrect when coupling 
vanishes after passing a coupling region and more PESs 
are occupied. In the MSDR this problem is often less 
significant, because in many cases after passing the cou- 
pling region fidelity does not decay anymore. The FSSH 
dynamics also suffers from several problems besides the 
inaccuracies caused by the classical description of nu- 
clei, such as the problem of "excessive coherence" related 
(similarly as for the LMFD) to the fact that all matrix 
elements of the density matrix attached to a trajectory 
evolve on the same PES. In many cases, this can be al- 
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leviated by applying the "decoherence" correction<i^— 
Including a similar correction into the MSDR method 
should be straightforward and will be explored in future 
work. Another possibility how to overcome some defi- 
ciencies of the LMFD or FSSH dynamics would be to use 
the MSDR together with one of the propagation meth- 



ods attempting to solve the mixed quantum-classical Li- 
ouville equation directly. 
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